home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / VideoToolboxSources / DrawGrating.c < prev    next >
Text File  |  1995-10-09  |  5KB  |  172 lines
  1. /*
  2. DrawGrating.c
  3.  
  4. Quickly draws a grating (cos, sin, square, or missing fundamental) with a
  5. gaussian envelope in the current port, bounded by the supplied rect. The routine
  6. is accurate and takes advantage of even and odd symmetries for speed. It is not
  7. entirely general, e.g. orientation is restricted to vertical and horizontal.
  8.  
  9. Squarewaves have an undefined value at zero crossings, e.g. at the origin, so I've
  10. shifted all the functions by a half pixel up and to the left.
  11.  
  12. WARNING: To compile this routine you'll need the header files nr.h and nrutil.h,
  13. which are part of the Numerical Recipes in C, sold by Cambridge University Press
  14. for about $100. (See Advice document.) Furthermore, within nr.h, you should make
  15. sure that ANSI prototypes are enabled, since they provide the strictest and
  16. safest prototype information. Similarly, set your compiler to REQUIRE prototypes.
  17. The FLOAT type is my modification of the Numerical Recipes files, as explained in
  18. the "Numerical Recipes.doc" file in Notes. (Bosco Tjan reports that you need the 
  19. current version, 2, since version 1 omits the lvector() function.)
  20.  
  21. HISTORY:
  22. 11/2/94 dgp shrink by s->noiseCheckSize, since it will later be expanded by that amount.
  23. 3/26/95 dgp as requested by Manoj, there are now separate space constants along and across
  24.             bars.
  25. 6/16/95 dgp replace "noiseCheckSize" by "signalCheckSize".
  26. 6/17/95 dgp changed "signalBounds" in first line to "stimulusBounds"
  27. 6/17/95 dgp eliminated all reference to s->signalContrast, since that applies to the 
  28.             entire sequence of images.
  29. */
  30. //#include "TextInNoise.h"
  31. #include "nr.h"            /* prototypes of Numerical Recipes and definition of FLOAT */
  32. #include "nrutil.h"
  33. enum{kCosinewave=0,kSinewave,kSquarewave,kMissingFundamental};// gratingType
  34. typedef struct{
  35.     FLOAT degreesPerPixel;
  36.     FLOAT spatialFrequency,spaceConstantAlongBars,spaceConstantAcrossBars;
  37.     Rect bounds;
  38.     Point origin;    // origin of the functions will be half a pixel up and to the left of this.
  39.     int gratingType;
  40.     int checkSize;
  41.     Boolean vertical;
  42.     int eAmplitude,eOffset;
  43. }GratingSpec;
  44. void DrawGrating(GratingSpec *g);
  45.  
  46. #undef MAX
  47. #undef MIN
  48. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  49. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  50.  
  51. void DrawGrating(GratingSpec *g)
  52. {
  53.     unsigned long *gratingPix;
  54.     FLOAT *fX,*fXEnvelope,*fY,*fTemp,w,e1,b,dc,a;
  55.     int iMax,width,fXSymmetry=0,fYSymmetry=0;
  56.     enum{kEven=1,kOdd};
  57.     Rect r;
  58.     int i,j;
  59.     int leftEdge,verticalOrigin;
  60.  
  61.     r=g->bounds;
  62.     OffsetRect(&r,-g->origin.h,-g->origin.v);
  63.     ShrinkRect(&r,g->checkSize,g->checkSize);    // Shrink
  64.     iMax=MAX(r.right,-r.left)-1;
  65.     iMax=MAX(iMax,MAX(r.bottom,-r.top)-1);
  66.     fX=vector(-iMax-1,iMax);
  67.     fXEnvelope=vector(-iMax-1,iMax);
  68.     fY=vector(-iMax-1,iMax);
  69.     width=r.right-r.left;
  70.     b=1/(g->spaceConstantAlongBars/g->degreesPerPixel);
  71.     b*=g->checkSize;    // Shrink
  72.     fYSymmetry=kEven;
  73.     for(i=0;i<=iMax;i++){
  74.         a=(i+0.5)*b;
  75.         fY[-i-1]=fY[i]=exp(-a*a);
  76.     }
  77.     if(g->spaceConstantAlongBars!=g->spaceConstantAcrossBars){
  78.         b=1/(g->spaceConstantAcrossBars/g->degreesPerPixel);
  79.         b*=g->checkSize;    // Shrink
  80.         for(i=0;i<=iMax;i++){
  81.             a=(i+0.5)*b;
  82.             fXEnvelope[-i-1]=fXEnvelope[i]=exp(-a*a);
  83.         }
  84.     }else{
  85.         for(i=0;i<=iMax;i++) fXEnvelope[-i-1]=fXEnvelope[i]=fY[i];
  86.     }        
  87.     assert(fYSymmetry==kEven);
  88.     w=2*PI*g->spatialFrequency*g->degreesPerPixel;
  89.     w*=g->checkSize;    // Shrink
  90.     e1=g->eAmplitude;
  91.     switch(g->gratingType){
  92.     case kCosinewave:
  93.         fXSymmetry=kEven;
  94.         for(i=0;i<=iMax;i++){
  95.             fX[-i-1]=fX[i]=e1*fXEnvelope[i]*cos((i+0.5)*w);
  96.         }
  97.         break;
  98.     case kSinewave:
  99.         fXSymmetry=kOdd;
  100.         for(i=0;i<=iMax;i++){
  101.             fX[i]=e1*fXEnvelope[i]*sin((i+0.5)*w);
  102.             fX[-i-1]=-fX[i];
  103.         }
  104.         break;
  105.     case kSquarewave:
  106.         fXSymmetry=kOdd;
  107.         for(i=0;i<=iMax;i++){
  108.             if(sin((i+0.5)*w)>=0)fX[i]=1;
  109.             else fX[i]=-1;
  110.             fX[i]*=e1*fXEnvelope[i];
  111.             fX[-i-1]=-fX[i];
  112.         }
  113.         break;
  114.     case kMissingFundamental:
  115.         fXSymmetry=kOdd;
  116.         for(i=0;i<=iMax;i++){
  117.             a=sin((i+0.5)*w);
  118.             fX[i]=(a>0)?1:-1;
  119.             fX[i]-=(4/PI)*a;
  120.             fX[i]*=e1*fXEnvelope[i];
  121.             fX[-i-1]=-fX[i];
  122.         }
  123.         break;
  124.     }
  125.     if(!g->vertical){
  126.         fTemp=fX;
  127.         i=fXSymmetry;
  128.         fX=fY;
  129.         fXSymmetry=fYSymmetry;
  130.         fY=fTemp;
  131.         fYSymmetry=i;
  132.     }
  133.     gratingPix=lvector(-iMax-1,iMax);
  134.     dc=g->eOffset+0.5;
  135.     leftEdge=r.left+g->origin.h/g->checkSize;
  136.     verticalOrigin=g->origin.v/g->checkSize;
  137.     for(j=iMax;j>=0;j--){
  138.         if(j>=r.bottom && -j-1<r.top)continue;
  139.         switch(fXSymmetry){
  140.         case kEven:
  141.             for(i=iMax;i>=0;i--)gratingPix[-i-1]=gratingPix[i]=dc+fY[j]*fX[i];
  142.             break;
  143.         case kOdd:
  144.             for(i=iMax;i>=0;i--){
  145.                 a=fY[j]*fX[i];
  146.                 gratingPix[i]=dc+a;
  147.                 gratingPix[-i-1]=dc-a;
  148.             }
  149.             break;
  150.         default:
  151.             assert(0);
  152.         }
  153.         if(j<r.bottom)SetPixelsQuickly(leftEdge,j+verticalOrigin,&gratingPix[r.left],width);
  154.         if(-j-1>=r.top){
  155.             switch(fYSymmetry){
  156.             case kEven:
  157.                 break;
  158.             case kOdd:
  159.                 for(i=iMax;i>=-iMax-1;i--)gratingPix[i]=dc-(gratingPix[i]-dc);
  160.                 break;
  161.             default:
  162.                 assert(0);
  163.             }
  164.             SetPixelsQuickly(leftEdge,-j-1+verticalOrigin,&gratingPix[r.left],width);
  165.         }
  166.     }
  167.     free_lvector(gratingPix,-iMax-1,iMax);
  168.     free_vector(fX,-iMax-1,iMax);
  169.     free_vector(fXEnvelope,-iMax-1,iMax);
  170.     free_vector(fY,-iMax-1,iMax);
  171. }
  172.